# =============================================================================
# Statistical Analysis Code for:
# "Chronic adaptive versus conventional DBS response patterns in Parkinson's
#  disease: A pilot randomized crossover trial"
#
# Sensitivity Analysis: Reduced ANCOVA and Wilcoxon signed-rank test
# Author: Jun Tanimura, MD, MSc.
# Created: 2025-11-13
# Last Updated: 2025-12-07
#
# Note:
# - Data preprocessing is shared with 02_ANCOVA_v3.1_FULL.R (main analysis)
# =============================================================================

library(lme4)
library(dplyr)
library(tidyr)
library(broom.mixed)
library(effectsize)

# Data loading and common objects
# Extract data preprocessing from main script

# Explicitly set Treatment effect sign: aDBS - cDBS by setting ref to cDBS
dat_long[[treat_var]] <- relevel(dat_long[[treat_var]], ref = "cDBS")

# Reduced ANCOVA (without crossover terms)

fit_reduced_ancova <- function(df, outcome, id_var, treat_var, base_var) {
  
  f_reduced <- as.formula(
    paste0(outcome, " ~ ", treat_var, " + ", base_var, " + (1|", id_var, ")")
  )
  
  m <- lmer(f_reduced, data = df, REML = TRUE)
  
  # Extract Treatment coefficient (aDBS - cDBS)
  summ <- broom.mixed::tidy(m, effects = "fixed", conf.int = TRUE)
  
  tr_name <- paste0(treat_var, "aDBS")  # Corresponding row name (check with print(summ) if needed)
  
  est_row <- summ %>% filter(term == tr_name)
  
  est  <- est_row$estimate
  lwr  <- est_row$conf.low
  upr  <- est_row$conf.high
  pval <- est_row$p.value
  
  # Cohen's d (for paired design). Simplified approach:
  # Calculate d from t-value and df, or use effectsize::t_to_d()
  t_val <- est_row$statistic
  df_res <- est_row$df
  d_obj <- effectsize::t_to_d(t = t_val, df_error = df_res, paired = TRUE)
  d     <- d_obj$d
  d_lwr <- d_obj$CI_low
  d_upr <- d_obj$CI_high
  
  tibble::tibble(
    outcome   = outcome,
    model     = "ANCOVA_reduced",
    estimate  = est,
    ci_lower  = lwr,
    ci_upper  = upr,
    d         = d,
    d_lower   = d_lwr,
    d_upper   = d_upr,
    p_value   = pval
  )
}

# Wilcoxon signed-rank test (simple paired comparison)

fit_wilcoxon <- function(df, outcome, id_var, treat_var) {
  
  wide <- df %>%
    select(all_of(c(id_var, treat_var, outcome))) %>%
    tidyr::pivot_wider(
      names_from = !!sym(treat_var),
      values_from = !!sym(outcome)
    )
  
  # Assume column names aDBS, cDBS (check with names(wide) if needed)
  wide <- wide %>%
    mutate(delta = aDBS - cDBS)
  
  wtest <- wilcox.test(wide$aDBS, wide$cDBS, paired = TRUE, exact = FALSE)
  
  # Mean difference and t-distribution based CI
  m_diff <- mean(wide$delta, na.rm = TRUE)
  se_diff <- sd(wide$delta, na.rm = TRUE) / sqrt(sum(!is.na(wide$delta)))
  df_res <- sum(!is.na(wide$delta)) - 1
  t_crit <- qt(0.975, df_res)
  ci_lwr <- m_diff - t_crit * se_diff
  ci_upr <- m_diff + t_crit * se_diff
  
  # Cohen's d (paired)
  d_obj <- effectsize::cohens_d(wide$aDBS, wide$cDBS, paired = TRUE, ci = 0.95)
  
  tibble::tibble(
    outcome   = outcome,
    model     = "Wilcoxon_paired",
    estimate  = m_diff,
    ci_lower  = ci_lwr,
    ci_upper  = ci_upr,
    d         = d_obj$Cohens_d,
    d_lower   = d_obj$CI_low,
    d_upper   = d_obj$CI_high,
    p_value   = wtest$p.value
  )
}

# Apply to all outcomes

sens_results <- purrr::map_dfr(
  outcomes,
  ~ fit_reduced_ancova(dat_long, .x, id_var, treat_var, base_var)
) %>%
  bind_rows(
    purrr::map_dfr(
      outcomes,
      ~ fit_wilcoxon(dat_long, .x, id_var, treat_var)
    )
  )

# Save results (for supplement)
readr::write_csv(sens_results, "results_sensitivity_analysis.csv")
saveRDS(sens_results, "results_sensitivity_analysis.rds")